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I. Work Accomplished 

A. Impact on Mars 

One major result of studying impact cratering on Mars is given in: 

1. Stewart, S. T., J. D. O’Keefe, and T. J. Ahrens, Formation of fluidized ejecta blankets on Mars 
by shock-melting subsurface ice. Science, submitted, 2002. 

Many Martian craters are surrounded by ejecta blankets which appear to have been fluidized 
[Barlow, 1988; Mouginis-Mark, 1978] forming lobate and layered deposits terminated by one or more 
continuous distal scarps, or ramparts. One of the first hypotheses for the formation of so-called rampart 
ejecta features was shock-melting of subsurface ice, entrainment of liquid water into the ejecta blanket, 
and subsequent fluidized flow [Carr et al., 1977; Mouginis-Mark, 1979]. Our work quantifies this 
concept. Rampart ejecta found on all but the youngest volcanic and polar regions, and the different 
rampart ejecta morphologies [Barlow et al., 2000] are correlated with crater size and terrain [Barlow 
and Bradley, 1990; Mouginis-Mark, 1979]. In addition, the minimum diameter of craters with rampart 
features decreases with increasing latitude [Kuz'min et al., 1988] indicating that ice laden crust resides 
closer to the surface as one goes poleward on Mars (Fig. 1). 



Fig. 1. Geographic 
distribution of the onset 
of rampart (fluidized 
ejecta) craters on Mars. 
Onset diameters are 
greatest near the 
equator, indicating the 
depth (~ 1 km) to 

interstitial ice in the 
crust is greatest in the 
tropics (after Squyres et 
al. [1992]) 


Our ability to interpret the cratering record to learn about water on Mars has been limited by our 
understanding of the physical processes which fluidize ejecta blankets and the competing hypothesis that 
rampart morphologies may be produced by interactions between the ejecta curtain and the atmosphere 
[Barnouin-Jha and Schultz, 1998; Schultz, 1992; Schultz and Gault, 1979]. 

In Paper 1, we report on the results from numerical simulations of impact cratering onto an ice- 
rich Martian regolith, including the projectile and ejecta flow through the atmosphere, and the 
implications for formation of fluidized ejecta and global crustal water content. 

Recently, we have experimentally demonstrated from our measured shock properties of H 2 0 ice 
[Stewart and Ahrens, 2002] that shock-induced melting begins at shock pressures of only 0.6 GPa 
under terrestrial conditions, more than an order of magnitude lower than the 6.2 to 10 GPa previously 
inferred [Ahrens and O'Keefe, 1985; Kieffer and Simonds, 1980]. Ice in the Martian crust, at 
temperatures between 150 and 273 K, will begin to melt after experiencing shock pressures between 2.2 
and 0.6 GPa, respectively, and will completely melt upon release from shock compression to pressures 
above 6 and 3.7 GPa. 

We conducted long duration simulations of crater formation by asteroid-derived impactors onto 
an ice-filled Martian regolith, starting with the projectile several atmospheric scale heights above the 
surface (30 to 60 km) and continuing through the formation of the ejecta blanket and collapse of the 
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impact vapor plume. Using the finite-difference, Eulerian-Lagrangian hydrocode CTH [McGlaun et ah, 
1990], normal impacts were simulated using two-dimensional cylindrical symmetry. The spatial 
resolution of the Eulerian mesh was constrained to 20-30 cells across the projectile radius, and the width 
and depth of the simulated crust were about an order of magnitude larger than the final crater size to 
minimize edge effects. Several hundred massless Lagrangian tracer particles were distributed through the 
crust to track the shock pressure decay profile with depth, ejection angles and velocities, and the 
excavated zone. 

The model crust and atmosphere were initialized as isothermal layers in gravitational equilibrium 
with a surface pressure of 7 mbar. The initial temperatures of the atmosphere and crust were 1 70 and 200 
K, respectively, which is representative of Martian mid-altitudes and mid-latitudes. We considered a 
range of ground ice content using an exponential decay profile to simulate a theoretical regolith with ice- 
filled cracks and pore spaces each computation cell representing the crust was initialized with the 
specified ice and rock volume fractions under lithostatic pressure, and the explicit equation of state for 
each material is utilized throughout the computation. By assuming equal shock pressures in the ice and 
silicate fraction, upon partial melting ice marked strength reduction follows as we volume averaged yield 
strength in the mixture. The atmosphere is considered a pure C0 2 ideal gas. Asteroidal impactors and the 
Martian crust are probably basaltic in composition. However, a shock equation of state for dunite 
ANEOS [Thompson and Lauson, 1972] was used as an approximation. In mixed cells, equations of 
state for each component, the C0 2 atmosphere, silicate projectile and crust, and H 2 0 ice, are used to 
calculate shock pressures and temperatures. The unmelted dynamic strength of the Martian crust is 
constrained by the observed transition diameter, D T , from simple to complex craters on Mars, ~7 km 
[Garvin et ah, 2002], Our previous cratering studies (Eq. 28, O’Keefe and Ahrens [1993]) define an 
effective yield strength of 10 7 Pa [O'Keefe and Ahrens, 1999; O'Keefe et ah, 2001] for Mars. We 
varied the projectile diameter (100-2000 m) to simulate simple, transition and complex type crater 
formation at typical asteroidal impact velocities, 10 km/s. 

The shock wave pressure decays exponentially with distance from the impact point [Melosh, 
1989], Using the pressure decay profile, we calculate the radial extent of impact-induced melting of 
ground ice. The flow field of the Lagrangian tracer particles determines the zone of excavation used to 
calculate the fraction of ice melted within the ejecta blanket. 

The results of the simulations are summarized in Fig. 2. The transient crater diameter corresponds 
to the hemispheric diameter at the maximum depth of penetration, and the final crater diameter is given 
by the location of the overturned ejecta flap. The crater diameters are in very good agreement with the 
crater scaling laws developed for competant rock surfaces. The depths for shock-induced partial and 
complete melting were determined from the shock-pressure decay profiles. Fig. 3. For rocky impactors at 
10 km/s, the peak shock pressures are just above 100 GPa, and very little rock or ice is vaporized. The 
bulk of material excavated into the continuous ejecta blanket originates from the zones of complete and 
partial melting. Fig. 4. For a near-surface ground ice content of 20% (by volume), we find that the 
hemispherical volume where subsurface ice is partially melted by the impact shock has a radius of about 
14 projectile radii (r p ), and ice is completely melted within about 7 r p . Partial melting of ground ice occurs 
throughout the excavation zone. The extent of melting will slightly increase as the composition 
approaches pure rock (Fig. 3), although the melt zone is expected to be slightly smaller for the more 
common oblique impact events [Pierazzo and Melosh, 2000]. Also, shock-induced melting may be 
more efficient if the ice has significant salt content. The presence of 0-20% (by volume) ground ice 
increases the diameter of the final crater by only a fraction of a projectile radius (Fig. 4A and B). 

The material ejected near the crater rim originates from both the completely and partially melted 
zones (Fig. 4), and more than half the ice which is excavated into the continuous ejecta blanket is melted. 
For 10 km/s impacts, the ejection velocities decrease from almost 3 km/s near the impact point to 10's m/s 
near the crater rim. Thus, the solid ejecta curtain following approximately ballistic trajectories will impact 
the surrounding terrain with horizontal flow velocities in the range of 10's- 100's of m/s, within the range 
of large terrestrial landslides. 
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Fig. 2. Transient ( X ) and final ( □ ) crater diameter as a function of projectile diameter for typical asteroid-derived 
impactors on Mars. Diameter of hemispherical zone of partial melting of ground ice (dashed line) is comparable to 
final crater diameter. Size of zone of complete melting (dotted line) is about half the transient crater cavity. Solid 
lines correspond to crater scaling functions in competent rock [ Melosh , 1989]. 



Fig. 3. Shock pressure decay profiles for different Martian crust composition. Horizontal lines indicate critical 
shock pressures for partial and complete melting and complete vaporization. 
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Fig. 4. Crater excavation zones for 
500 m diameter projectile impacting 
at 10 km/s into (A) pure rock and 
(B) a 20% (volume) ice-rock 
mixture. Trajectory histories shown 
for Lagrangian tracer particles ( X ) 
from within the excavated crater 
cavity. Ejecta trajectories shown for 
computation with no atmosphere at 
20 sec. 



Fig. 5. Time sequence for 500 m 
diameter projectile impact onto 
regolith with 20% (volume) near- 
surface ice, including Martian 
atmosphere. Each panel shows 
ice content (left, shading), 
temperature contours (left), flow 
field velocities (right), and major 
material composition (solid 
colors). 


Fig. 5 displays a time series from a simulation with 500 m-diameter projectile, forming a 7.5-km 
diameter crater, showing the atmospheric shock (A), ejecta curtain, and interaction with the vapor plume, 
and collapse of the impact vapor plume. For simple craters, with D < 7 km, the water vapor plume is 
contained by the mass of the atmosphere and does not spread far beyond the crater rim. For large craters, 
D>10's km, a portion of the water vapor mass may be dispersed into the upper atmosphere, transiently 
affecting the regional and perhaps global climate, and remaining mass of the plume collapses down over 
the crater and ejecta blanket as transient rain or sleet. For transition size craters, D~7 km, the vapor plume 
appears to interact with the solid ejecta curtain during the formation of the continuous ejecta blanket. 
Figs. 6-8 summarize the simulations for a 0. 1 , 0.5, and 2 km diameter silicate impactors. 

Models for the formation of rampart morphologies by fluidized flow [Ivanov, 1996] suggest a 
conservative liquid water volume fraction, corresponding to a near-surface ice volume fraction of about 
20-40%. Using the observed rampart ejecta craters [Barlow and Bradley, 1990], corrected for 
observational bias using the Martian crater production function [Hartmann, 1999]. Summing over the 
entire Martian rampart crater size range, D=l-32 km, the volume of ice implied by the observed Martian 
rampart crater number distribution is 20% which is equivalent to a global layer of water about 0.6 m deep. 
Based on the crater production function, this value corresponds to the amount of ice excavated by impact 
craters over the past 3 Ga. Since rampart craters stochastically sample a fraction of the surface, limited to 
excavation depths within the upper 2 km for the observed size range, the actual regolith ice content must 
be considerably larger. Only 0.5 % of the Martian surface area was sampled by the observed rampart 
craters. Correcting for the sampling area, the implied total ice content, within 2 km of the surface, is 
equivalent to a global layer of water about 120 m thick. 


4 



B. Impact Damage and Porosity Cratering Studies 

Our major results in this area are in: 

2. O’Keefe, J. D., S. T. Stewart, M. E. Lainhart, and T. J. Ahrens, Damage and rock-volatile 
mixture effects in impact crater formation, Int. J. Impact Engng, 26, 543-553, 2001. 

3. O’Keefe, J. D. and T. J. Ahrens, Hydrodynamic instabilities and impact cratering, submitted to 
Nature [2002]. 

Our goal in Paper 2 was to determine what strength model(s) reproduce the faults and complex features 
found in large scale gravity driven craters. Collapse features found in large scale craters require that the 
rock strength weaken as a result of the shock processing of rock and the later cratering shear flows (e.g. 
Melosh and Ivanov [1999]). In addition to the presence of molten silicate in the intensely shocked region, 
the presence of water, either ambient, or the result of shock melting of ice weakens rock. There are several 
other mechanisms for the reduction of strength in geologic materials including dynamic tensile and shear 
induced fracturing (e.g. Melosh and Ivanov [1999]). Fracturing is a mechanism for large reductions in 
strength. We found that by incorporating damage into the models that we could in a single integrated 
impact calculation, starting in the atmosphere produce final crater profiles having the major features found 
in the field measurements [. Morgan et al., 1997] (central uplifts, inner ring, terracing and faulting. This 
was accomplished with undamaged surface strengths (~0.1GPa) and in depth strengths (~1.0 GPa). 

We modeled the damage in geologic materials using a phenomenological approach, which 
coupled the Johnson-Cook damage model [Johnson and Holmquist, 1994]. with the CTH code 
geologic strength model. The objective here was not to determine the distribution of fragment sizes, but 
rather to determine the effect of brecciated and comminuted material on the crater evolution, fault 
production, ejecta and final crater moiphology. There are a number of definitions of damage which are 
discussed by Kachanov [1986], Damage (D) here is defined as the ratio of the integrated plastic strain 
experienced by a given material particle relative to the integrated plastic strain required at failure (loss of 
strength). The time derivative of damage is given by 

dD/dt = ( de/dt ) /e f ( 1 ) 

where de /dt is the plastic strain rate. The integrated plastic strain at failure, e f is defined 
by 

e f = ( e fo + Efduc exp( - £„, P/ Y 0 )*( 1 + £ f d£/ dt In (de/dt) )*( 1 + £„ (T-T r )/(T m -T r ) ) (2) 

and e ffi , e fduc , £(p , e fd£/dt # and Err . are material parameters, P and T are the local pressure and temperature 
an d T m and T r are the melt temperature and the thermal softening temperature, Y 0 is the undamaged Mohr- 
Coulomb yield strength at zero pressure, £ ro is a measure of the magnitude of the strain at failure at 
ambient temperature and zero pressures. £(p is a measure of the rate of increase of ductility with pressure 
and hence the transition from brittle to ductile behavior [Scholz, 1990]]. £ fdE/dt is a measure of the rate of 
change in £ f due to strain rate effects, and is important for small scale impacts, and is not as important for 
large scale late stage shear flows. % is a measure of the increase in ductility due to thermal effects. The 
above parameters have been determined for a number of metals , concrete and some ceramics. Rocks 
have had little work. Recently, more laboratory damage data has become available on geologic materials 
on account of efforts of researchers at LLNL [Lomov et al, 2001] and LANL [Swift et al., 2000].. 

The damaged strength at zero pressure is a function of thermal softening, and damage and is 
given by 

Y 0d = Y dT (l-D) (3) 

where Y dT is the thermally and density degraded strength at zero pressure. The damage function is unity 
at failure and zero when there is no degradation. This reduction in strength was incorporated into the 
geologic strength model to account for pressure effects (Mohr-Coulomb ) and is given by 

Y s = Y vm + (Y 0d -Y,J exp { -dY s /dP* P /( Y^-Yo,) } (4) 

where Y^ is the limiting strength at high pressure (von Mises magnitude). The internal angle of friction 
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Mars Simulation: Impact onto ice-rich surface with atmosphere 
Formation of simple crater. Simulation duration 60 sec. 

Rocky projectile: R=50 m V=10 km's Zo=33 km 
Target contains 20%vol ice at surface decreasing 
with depth by exponential constant of 3 km. 

Isothermal crust 200 K. effective strength 10 MPa. 

Pure CO 2 ideal gas almosphere: 
isothermal 1 70 K. surface pressure 7 mbar. 
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Fig. 6. Particle velocity, material field and ice contents for 50 m radius silicate 
impactor intersecting with CO 2 atmosphere, 20% icy surface of Mars @ 10 km/sec. 
(a) Time = 3.25 sec projectile 1-km above surface, (b) Time = 26.5 sec. water 
vapor plume has expanded over transient cavity; 1.8 km in diameter, (c) Time = 
27.5 sec water vapor plume starting to collapse, (d) Time — 60.0 sec vapor plume 
has collapsed and Mars’ atmosphere is reestablished above 1.8 km diameter crater 
that has not been modified. 
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Mars Simulation: Impact onto ice-rich surface with atmosphere 

Formation of transition size crater. Simulation duration 180 sec. 

Rocky projectile: R=250 m V=10 km/s Zo=61 km 
Target contains 20% V ol ice at surface decreasing 
with depth by exponential constant of 3 km. 

Isothermal crust 200 K, effective strength 10 MPa. 

Pure CO2 ideal gas atmosphere: 
isothermal 170 K. surface pressure 7 mbar. 
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Fig. 7. Particle velocity and material field and ice contents for 250 m radius impactor impacting Mars, (a) Time = 5.9 
sec, projectiles is 2 km above surface, (b) Time = 52.5 sec. Transient crater, 6 km iin diameter formed and vapor 
plume has risen to ~50 km altitude, (c) Time “113 sec. Transient cavity is slightly relaxing as water plume 
collapses, (d) Time — 1 80 sec. Crater slightly collapses as CO 2 atmosphere is established above crater. 
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Mars Simulation: Impact onto ice-rich surface with atmosphere 

Formation of peak-ring crater. Duration 300 sec. 

Rocky projectile: R=1 km V=10 km/s Zq= 70 km 
Target contains 20% VO | ice at surface decreasing 
with depth by exponential constant of 3 km. 

Isothermal crust 200 K, effective strength 10 MPa. 

Pure COg ideal gas atmosphere: 
isothermal 170 K, surface pressure 7 mbar. 
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Fig. 8. Particle velocity and material 
field and ice contents for 1 km radius 
impactor, impacting Mars, (a) Time - 7 
sec. Projectile is just striking planet, (b) 
Time - 29 sec. Transient cavity formed, 
18 km in diameter and 100 km vapor 
plume driven into atmosphere, (c) Time 
- 77 sec. Transient cavity is collapsing, 
(d) Time - 70 sec. Central peak extends 
5 km upwards, (e) Time - 258 sec. 
Central peak has collapsed leaving 
ringed crater structure ~ 40 km in 
diameter. Local atmosphere remains as 
vaporized water gas from regolith. 






also degrades with damage and is given by 

, dY/dp = dY/dP] Dam - (1-D){ dY/dPWdY/dP] Dam } (5) 

where dY/dP] m Dam and dY/dP] Dam are the imdamaged and damaged measures of internal friction. The 
above yield surface is plotted in Fig 9. In the case of no damage, and at high pressures, the yield strength is 
the von Mises limit. At zero pressure , the yield strength is given by Y 0(J (eqn 4), and under extensive 
conditions, failure (zero yield strength ) occurs at the tensile failure limit, -T f . In the case of complete 
damage, the yield strength is zero at zero pressure and increases with confining /overburden pressure, at a 
rate depending upon dY/dP ( eqn. 5 ). The yield surface limit at high pressures is given by the damaged von 
Mises limit, Y vmd . The shape of the damaged yield surface depends upon the degree of damage to the 
internal angle of friction, dY/dP ( eqn. 5 ). Examples of damage zone for an impact representative of the 
Chicxulub event and the effect of overburden pressure on the local strength during that impact is shown in 
Fig 10. 

We accomplished our goal of determining whether or not the damage approach to was capable of 
reproducing faults and complex features found in large scale craters with plausible assumptions for the 
magnitudes of the material parameters. In general, we only found faulting when damage was modeled. 
Specifically, we performed parametric variations on the magnitudes of some of the key material constants 
to determine if that resulted in plausible values for e f and dY/dP and to determine if these would result 
features such as in inner rings, central peaks and faults. In the expression for £ f , we varied £«, and set the 
other parameters to zero. This bounds the magnitude of £ t ~ . since the other terms in general increase the 
magnitude of e f .We found that for the Chicxulub impact case, that the parameter ranges for central uplift 
and inner ring formation are when the in depth strength (Y vm ) < l.OGPa, the internal friction angle <45 
deg. and the £ f is <0.15. We found for the conditions for fault scarps and slumping/terracing to form are 
when the impactor diameter is greater than 15 km , the damaged value of the internal angle of friction < 30 
deg. and £f ~ < 0.08. These integrated strain to failure values for both collapse and fault features are within 
the broad range of laboratory measurements. Strain at fault measurements of rock specimens under 
confining pressure can vary from 0.03 to 0.25 (0.2GPa confining pressure). Rocks within a large scale 
impact environment can be expected to have values >0.1, because of shock heating and environmental 
thermal and pressure distributions [Brace et al. , 1 966] . 

Dilatancy has long been recognized as an important aspect of rock response [Brace et al., 
1966], it is produced by the creation of cracks from inhomogeneous stresses during loading and at low 
pressures it results in a porous, decreased density, rock. (Fig.l 1). We have developed a first order model 
for bulking. The model is a relationship between the degree of damage, pressure, and density. This is 
shown in Figure 11. For imdamaged rock, the density increases inversely with the magnitude of the bulk 
modulus. For completely damaged rock, which can be characterized as gravel, the density is 20 to 40% 
less than undamaged or competent rock. We assume that the gravel density is the limiting density of 
damaged rock at zero pressure. The bulk moduli has been measured for many gravels and its magnitude is 
much less than the competent rock The density of the damaged rock(gravel) also increases with the 
inverse of the gravel bulk modulus. With increasing pressure, the cracks close and the material locks 
together at a pressure determined to be the intersection of the bounding density curves. This point is 
shown in Fig 11. We also assume that the damaged density at zero pressure varies at linearly with 
damage and the magnitude of the calculated bulk modulus requires the density- pressure curve to intersect 
at the bounding intersection point, (see Fig 1 1). An example of the effect of bulking on crater profile is 
shown in Fig. 12. The decrease in the crater depth correlates directly with the decrease in rock density 
due to fracture shown in Fig. 11. 

In our recent survey (Paper 3) of hydrodynamic instabilities in various planetary hypervelocity 
impact scenarios, we have explored a very wide range of impacts of both similar and strongly dissimilar 
projectile/target density ratios, DR = 5/p where 5 is projectile density and p target density and VN is the 
vaporization number U /E v . What we found is that when the specific kinetic projectile energy far exceeds 
the target vaporization energy target instabilities always form (Fig. 13). The regime previously explored 
by Holsapple and Schmidt [1982] in which the projectile spreads and smoothly lines the transient crater 
cavity occurs only in the case of VN < 10 and DR < 3. For greater values of VN and DR, projectile or 

target/RT (Rayleigh-Taylor) or K/H (Kelvin-Helmholtz) instabilities grow faster than the slowing of the 
projectile by the target material. Thus these instabilities control transient crater growth 
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Fig. 9. 


Fig. 10 


Fig. 9. Diagram of the yield surface as a function of normalized pressure. Top curve is case of no damage (D=0.0) 
and bottom case is complete dqamage (D=l .0). 

Fig. 10. An example of the effect of damage on the yield strength distribution. Note the strength in the damaged 
cavity due to the overburden pressure. The magnitude of this yield strength depends upon dY/dP. If this is zero, then 
there is no strength in the damaged cavity and the material acts like a liquid, and would not exhibit features such as 
inner rings and central peaks. The time is 40 s and is at the time of maximum penetration. Impact is representative 
of Chicxulub event, velocity = 20 km/s, diameter = 12 km. 


and the final crater profile. Moreover, or VN < 10 and DR > 3, the ejecta mass and velocity decrease 
with DR, with no ejecta for DR >10. 

An example of the changes in the key final crater measures with increasing DR is shown in Fig. 14. 
The changes in the density ratio was accomplished by varying the density of the projectile or by varying the 
porosity of the target, which in both cases was gravity with a strength of 0. 1 GPa. The present calculations, 
predict that for a solid granite projectile impacting a solid granite target a normalized volume, V/a 3 , of 638. 
Holsapple’s formula predict a volume of 264. Moreover, the present calculations predict a much faster 
increase in normalized volume with increasing values of DR (Fig. 14) than is predicting by similarity 
scaling theory and also discontinuous changes in final crater measures (Fig. 14). 


Determined from 
measurements 



Fig. 11. Loading paths for dilatancy. Model is constrained 
by properties of consolidated rock at zero damage and 
granular material at complete (D=1.0) damage. Material 
would initially load along its D = 0.0 line until completely 
damaged, and if completely damaged, unload along the 
release adiabat (D=1.0) line. 


H/a 



Fig. 12 Final crater profiles for cases showing the 
effect of dilatancy [ O'Keefe and Ahrens , 2000]. 
Note decreased crater depth with bulking. The 
impact dimensionless parameters are Y m /pU 2 = 
2x1 O' 4 and ga/U 2 =3.4 x 10' 5 , where U = impact 
velocity, p = target density, a = radius. 
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